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INTRODUCTION 


This manual comprises the final report on NASA Grant NGR 06-002-191. 
It describes a finite element computer code for the analysis of mantle 
convection. The coupled equations for creeping viscous flow and heat 
transfer can be solved for either a transient analysis or steady-state 
analysis. For transient analyses, either a control volume or a control 
mass approach can be used. Non-Newtonian fluids with viscosities which 
have thermal and spacial dependencies can bo easily incorporated. All 
material parameters may be written as function statements by tne user 
or simply specified as constants. A wide range of boundary conditions, 
both for the thermal analysis and the viscous flow analysis can be 
specified. For steady-state analyses, elastic strain rates can be 
included. « Although this manual was specifically written for users 
interested in mantle convection, the code is equally well suited for 
analyses in a number of other areas including metal forming, glacial 
flows, and creep of rock and soil. 
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GOVERNING EQUATIONS 

In this section the governing equations for the coupled thermal 
mechanical behavior of creeping viscous flows are presented, They are 
provided primarily for definition, therefore only a minimum of discussion 
is given. 


Creeping Elasto Viscoplastic Flow 
Equilibrium: 


°ij,j + s X i - 0 


CD 


Incompressibility : 


Velocity, strain-rate and vorticxty: 


e. . = — ?u. . + u. .) 
13 2 1,3 3,1' 



C 3 ) 

w. . = ^ (u. . - u. .) 
13 2 3,1 1,3 



C 4 ) 

Constitutive : 




c . , — o J , + a]. V 

13 2 \i 13 2 G 13 



C 5 ) 

„ So. . 8a! . 

0 . + u __±1 . 
13 3 t u k 

a! u> . 
ip P 3 

- a! a . 
3 P Pi 

C6) 

13 13 * 13 



C 7 ) 

P 58 - T 0 . . 
r 3 11 



C8) 

0 . P 

s . . = e . . + e , . 

13 13 3-3 



C 9 ) 

Boundary Conditions: 




a. .v. = T. = C, . (u. - u.) 
3-3 3 1 13 3 3 

on S 

u 


CIO) 

a. .v. = T. = f. 
13. 3 1 1 

on S 

0 


CH) 
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All torms arc defined in the tablo presented at the end of this section. 
It should bo noticed that the equations for pure v'seous flow are 
obtained by lotting G « . 

In the abovo equations, the viscosity can be a function of 
tomporature and the invariants of the stress tensor. Often it is 
assumed that tho second invariant of the stress deviator is the only 
stress invariant of importance. The following equations define this 
invariant (the effective stress) and its relationship to the effective 
strain rate. 


0 


off 


ij 


a 'u 


( 12 ) 


e eff 




(13) 


Tho relationship between tho viscosity and the effective stress 
and the effective plastic strain rate can bo found as follows. 


a!. * 2\i e?. 
ij ij 

(14) 

a ! . 0 ! . = 2\i c? . 2 y £? . 
13 13 13 13 

(15) 

2 2 , 2 3 2 

3 a eff 4m 2 e eff 

(16) 

0 eff 

M * 3e p 
^eff 

(17) 


It should be noted that for a state of uniaxial stress and incompressible 
plastic flow the effective stress and effective plastic strain rate 
reduce to the simplified form: 


0 eff " a axial 


( 18 ) 



'axial 


(19) 



This convenient relationship allows any one dimensional constitutive 
theory to be interpreted as a relationship between the effective stress 
and effective plastic strain rato. 


lloat Transfer 


Conductive, convective hoat transfer: 


n r r*£ + „ ILl 

pC p Lst i 3x.J 


1 n 3 [ 

k ill 

J [ 

, k 


+ Q 


Boundary conditions: 


K V - 11 Mr] \ - ‘i* 


or 


k iiv. 
3x. i 
1 


where 
q* = q* \ 
or 


= q ‘ 


on 


S 

q 


q = q / 

and 

q = h(9 - ej on s Q 
and 

v. = unit outward normal vector to S 
i 

Viscous heating: 


( 20 ) 


( 21 ) 


( 22 ) 


( 23 ) 


( 24 ) 
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Tablo of Terras and SI Units 


Symbol 

Definition 

SI Units 

C ij 

Surfaco viscosity tensor 

N-S/m 3 

c p 

Specific boat 

J/kg*k 

It 

Coefficient of convective surfaco heat 
transfer 

J/m 2 *S*k 

k 

Coefficient of thermal conductivity 

J/m*S'k 

P 

Prossurc * nogativo raoan normal stress 

N/m 2 

q 

Heat flow out of volume through surfaco 
duo to conduction 

J/m 2 *S 

q 

Known q on surfaco 

J/m 2, S 

q* 

Meat flow out of volume through surface 
due to conduction and mass transport 

J/m 2 *S 

q* 

Known q* on 

J/m 3 *S 

Q 

Heat source per unit volume 

s 

q 

Segment of surface whore q on q* 
is known 


s 

u 

Segment of surface whore velocity 
is known 


S 6 

Segment of surface where temperature 
is known 


S 

0 

Segment of surface where surface traction 
is known 


T i 

Surface traction 

N/m 2 

h 

Known T. on S„ 
1 o 

N/m 2 

u. 

1 

Velocity vector 

m/s 

“i 

Known u. on S ti 
1 u 

m/s 

E ij 

Total strain rate tensor 

1/s 

0 

s ij 

Elastic strain rate tensor 

1/s 

G P. 

XJ 

Plastic strain rate tensor 

1/s 
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Symbol Definition SI Units 

c 0 ££ Effective strain rate 1/s 

e eff Effective plastic strain rate 1/s 

0 Temperature 0 K or "C 

0^ Ambient temperature outside of volumo °K or °C 

2 

p Viscosity N*S/m 

v. Outward unit normal vector to surface 1 

P Density Kg/m 3 

2 

cJj.j Stress tensor N/m 

2 

Devlatoric stress tensor N/m 

2 

C eff Effectivo stress N/m 

Vorticity 1/s 


PjtEGEDJNG PAQ| 


BLANK NOT fILMkD 


VARIATIONAL PRINCIPLES 
and 

FINITE ELEMENT EQUATIONS 


The following section is included to provide a short summary of 
the finito element equations usod in MANTLE. For a tiK're complete 
development of these equations the reader should consult the references 
given in the Bibliography. 


Creeping Elasto Viscoplastic Flow 

The finite element equations develop from the following variational 
statement: 


6J ■ j j^ E ijdV 

V 


p6f. u dt- - j e a «p dv 

V V 


- f TjdUjdS - J pX i 6u i dV * 0 (26) 

S V 


where 


T i ■ Vi ■ C ij (u j - V on S u 


T. » T. 
i x 


on S 


(27) 


The above expression can bo derived from the governing equations 
via Galerkin's method or written directly as a statement of virtual 
work with the constraint of incompressibility incorporated through a 
Lagrange multiplier, p. 

Substitution of the constitutive equation presented in the previous 
section gives us 

6J = | 2p e.-fie.-dV - j £ ajJ fic-.dV - J pfi^.dV 

V V V 

- j e>, i 5p dV - j T^u-dS - | pX^U-dV = 0 

V S V 


( 28 ) 
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Wo now introduce the following natrix notations and finite element 


approximations within an element oi 

u A - tu(x,y)} - [N U ]{U> 0 (29) 

p » p(x,y) - [N p ]{p} o (30) 

G i;j “ {e(x,y?} « [N«](U} c (31) 

°ij * ta 'C?£»y)} 

oj[j « {cf' V (x,y)} 

{^(x^)} * [D]{e(x f y)> - IXKo ’ 7 (x»y) / (32) 

[D) b viscosity matrix 

[R] <* rolav’.ion matrix 

e ijL - [h]{e> (33) 

T i = [C]{u} c - [C]{u) e on S u (34) 

T. * CT> on S ff (35) 

X i * W 


Using this notation, tho contribution to 5J from element o can 
be written as 

M 0 ■ tu) T [K] e tu> e - {Su^[G]^(p } 0 - {«p)^[G]{u ) 0 

- tdu) T o {F) c - ♦ t«u)J[K c ] e tu) s (36) 


where 


[K ] 0 = 


[Ny] [D) [N f ]dV 


V(e) 




V(e) 


[N p ]‘[h] [N;]dV 


( 37 ) 

( 38 ) 
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{P} 


0 


* [ 
S>) 


{F C } 


e 


[K c ] 


V(o) 

S^Ce) 


[N u ] T tT>dS + j [N u ] T [C]{u} 0 dV 
S u (o) 

+ | P [N u ] T {X} dV 
V(o) 

[N^3 T [R)fa' V } dV 
[N u 1 T [(J] dv 


( 39 ) 

( 40 ) 

( 41 ) 


Summation of the clement matrices and enforcement of <SJ * 0 
for ail (6u) and (fip) provides the governing finite element 
equation 


K 

i t1 

1 G 
1 


u 


F 


P E 

G 

To" 


1 p 


0 

+ 

1 

0 


where [K] * l [K] e + l [Kj^. 


(42) 


Notb that as the coefficients in [C] become large, the above equation 
will force (u) = (u) on S . When this is the desired boundary 
condition on S , it can be specified directly rather than specifying 
a large C . . . 

1 j 

The code has also been written so that the penalty function 
approach can be used. This is accomplished using 


a. . 
xj 


* 2pe 


xj 


+ Xe, . 6 . . 

kk 


where X is specified as 10,000 * p . The final governing equation 
under this condition has the form 

[K]{u> = (F) + (F E ) . (43) 
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Hoat Transfer 

The finito dement equations for heat transfer are derived from the 

following variational statement f.und by Galorkin's principle, 

* 

<5J « 

4 

V 

- j q 5$ dS =0 (44) 

S 


34> . U_ 

3 *i 9X i 


Q 5$ + p C 


U 

3t 


1* * jf- «* 


dV 


q 


Substituting the finite element approximation for the temperature 
field and following the usual matrix summations, element by element, 
one obtains the governing equation 


+ [B]{*> = {Q} 


where 


[H] 


[N 1 ] T [R] [N’]dV + f [N)pC {U} T [N'] dV 
J P 

V V 


m 


[B3 



[N] T pC p [N] dV 
V 


{Q> 


[N] T Q dV - | [N) T q dS 

v s q 


4> a [N]{$> 


{V4>} 

q 


[N'K*> 



0 5) 


The heat flux across can be specified either as the total heat 

flux due both to convection and conduction, or the heat flux due to 
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conduction only. When is coincident with a streamline, there 
is no difference between the two specifications. However, if a control 
volume approach is used where mass is being transported across S , 
there is a significant difference between the two approaches. Specifi- 
cation of heat transfer by conduction only is tantamount to specifying 
the normal component of the gradient of the temporature. Because this 
is the more likely boundary condition to be known, MANTLE is written to 
accept either the temperature or heat transfer by conduction as known 
conditions. Note also that the coefficient of thermal conductivity 
can be specified for both the x-diroction and the y-direction. 



15 


PROGRAM MANTLE 

Introduction 

The namo given to the ovorall computer code is MANTLE. It is 
composed of four major subprograms (overlays on CDC computers) which 
are called by an executive program, DRIVE, The four programs aro, 
MESH1, a mesh generating program for mantle geometry; MESI12, a more 
general mesh generator; WAVE, a program which plans the solution pro- 
cedure and auxiliary storage of the matrix equations associated with 
the finite element analysis; and COUPLE, which is itself an executive 
program which calls programs TEMP and CREEP for thermal and plastic 
flow analysis, COUPLE also monitors the transient and nonlinear analy- 
sis and calls all output subroutines. Flow charts are provided in the 
following section. 

Each of the above programs will be discussed separately in terms 
of their input variables and whatever explanation is needed ' n define 
these variables. All input variables are echoed immediately after 
input with an output FORMAT identical to the input FORMAT. 


eregedibq pm* blan * not 
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Flow Charts 

Ti. following flow charts aro presented to clarify the use of 
control variables and convergence limits. They also indicate where in 
the program certain calculations are made, o.g., CPU time, strosses, 
or mosh adjustment during a quasi-Lagrangian analysis. However, they 
are far from being a complete description of the finite element code 
and in some instances imply a program flow different from that used in 
the actual code. For instance, the flow charts indicate the stiffness 
matrices are completely formed and then solved, rather than the 
piece-wise procedure actually used with the frontal method. These 
discrepencies do not distract from the usefulness of the flow charts, 
and in fact were allowed so as to emphasize the basic logic of the 
analysis. However, the reader should proceed with reasonable caution 
if he uses them to help him decipher the FORTRAN Code. 

woT 


t2 ss&fi& rtsa 
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back nubrttltutlon 
{DU> o [u]~ l [Dr 1 [l»]" 1 {pP} 


DELU o max | {du} f 
{U} - {u> + {DU} 


CALL 

BNDRYC 


CALL SECOND 
ITERC a ITERC + 1 
LCU = 0 


PRINT 

INCH, ITERC, DELU, 
DELP, SECOND 


"ITERC ;>ITMAXC 


'DELPHI 
DFCONV 
\ „ 


ITERC < 
INCLCU(2) 


INCLCU(2) = 
INCLCU(2) + INTLCU(2) 


LCU = 1 
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Tapo Usage 

A total of eight tapes aro designated in MANTLE. If all three 
major subprograms (MESII1 or MESII2, WAVE, and COUPLE) aro run as a unit, 
all eight tapes can bo scratch files. However, if several different 
runs arc to bo made with the same mosh, a significant savings in compu- 
ter timo will result if only COUPLE is called. This requires equating 
the tapes to permanent files. The following is an explanation of tho 
usago of each tape. 

TAPE1: ' This tape has tho output from WAVE related to the 
tapo segments for forming the stiffness matrix in 
CREEP. It must be a permanent file if this information 
is to bo used for later runs without calling WAVE, 

TAPE2: This tapo contains the segmented stiffness matrix for 
CREEP analysis. Because the matrix will bo in its 
LDU form, it can be saved for later analysis to 
avoid a new decomposition of tho stiffness matrix 
whenever it remains the same between runs, 

TAPE3: This tape is the counterpart of TAPE1 for TEMP. 

TAPE4: This tapo is tho counterpart of TAPE2 for TEMP, 

TAPES: All READ statements in MANTLE are from this tape. It 

should be oquatod to the INPUT file when data is read 
from punched cards. 

TAPE6: Normal output from MANTLE is placed on this tapo, 

that is, most WRITE statements are to this tape, 

It is normally equated to the OUTPUT file. 

TAPE7: Output which is useful for additional plotting and/or 

initialization for later runs in either time dependent 
problems or nonlinear analysis is read to this tape. 

If punch cards are used, it should be equated to the 
PUNCH file. Output from either MESH1 or MESH 2 is 
read onto this file by specifying I PUNCH = 7. 

TAPED: This tape contains all output from the MESH programs 

and WAVE not already on TAPE1 and TAPE3. It should 
be a permanent file if MESH and WAVE are not to be 
called for later runs. 
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Easy Input Roforenco 

The following listing of tho Input data is gi/cn ff' T ' tho user's 
quick reference. This listing is very abbr^Uted and primarily servos 
tho purpose of reminding an already familiar usor what variables are 
called for and in what order. 


PROGRAM DRIVE 


IFLOW(I> 


3110 


PROGRAM HESH1 


NPPE 

RlfRQrRMfRPI 

IPUNCH 

NBIVTHrNDIVR 
NPDCfXDCi YDCtTDC 
CHfCXfCYfTXfTYiTG 
NPDC r XDC r YBC r TBC 
CHfCXfCYfTXfTYfTn 
NIIMDC 

NPfNPBCrCOSXXPfXBCr YDCrTBC 
CHrCXfGYf TXrTYf TO 


no 

4E10.3 

4E10.3 

2110 

I 10r3E10*3 
6E10.3 
IlOf 3E10.3 
J10i3E10,3 

no 

2I10HE10t3 

4E10.3 


PROGRAM ME9H2 


NPPE no 

NUMLPS < 1 ) f NUMLPS ( 2 ) 2110 

XHINf XMAXrYMINfYMAX -1E10.3 

IPUNCH no 

NDIV(Ifl) rNDIV(If2) 2110 

JOIN( I r JrK) 4(I7rI3> 

XCQR ( I r J) r YCOR< I f J) 2E10.3 

NUMDC no 

J1 r LPBC t COBLDC > XLBC. YLBC t TOLBC 21 10 f 4E10 , 3 

CHDCrCXLBCf CYUBO»TXUDC»TYUBCf TOLBC 6E10.3 

NUMDC no 

Ilf NPBCrCOSXXPrXDCr YBG fTBC 2I10f3El0.3 

CHrCXf CYf TXiTYr TQ 6E10.3 


PROGRAM COUPLE 


NUMAT 

HATCI) 

INCPRt INCPUi INCPL 
INTPRf INTPU» INTPL 

ITV* MOP f TRANS r THETA r INTEMP r LAGEULf IRZ 
TIMER f HNI fDUMAXfETMAX 
ITMAXCfITMAXT 
VECTLf CTEMP 

INTLCU < 1 ) t INTLCU < 2 > r INTLTU { 1) t INTLTU< 2 ) rLCUrLTU 
DFCONV f DflCONV r DUCCNMr DTCONV 


no 

1015 

3110 

3110 

2.tl0r2E10.3r3I10 

E10.3rI10p2E10.3 

2110 


2E1U.3 

4110 

4E10.3 


NPTSf NSEC 

IlrXaRIUI>iYORD(I>fUX<I>fUYa)fUT<I> 
JBGN p JEND f INCR f TEMPO t UXO r UYO 


2110 

15r2E10.3f3ElB.10 
31 10f 3E10 • 3 
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Input Description 

PROGRAM DRIVE-- Input 

This is the main oxocutive program whose only purpose is to call 
the othor components of MANTLE and define the COMMON blocks which 
transfer data between those components. 


INPUT: 
/ — 


READ (5, 2) (IPLOW(I ) t I « l f 3) 
FORMAT CSI10) 


Those throe variables determine which of the four programs will 

be executed. The following code is employed; 

IF LOW Cl) B 1 calls MESH1 

■ 2 calls MESH 2 

IFL0WC2) - 1 calls WAVE 

IFL0WC3) = 1 


calls 


COUPLE 
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PROGRAM MBSH1-- Input 

Program MESH1 is a specialized mesh generator for mantle geometry. 
It is also through this program that the boundary conditions are speci- 
fied, although they can be altered lator in the analysis. Figure 1 
illustrates some typical meshes. The coordinate system used is 
rectangular rather than polar to accommodate the finite element approxi 
mation for pressure and its relationship to the derivatives of the 
velocity components. However, MESH1 accepts boundary conditions in 
terms of radial and tangential components. The mesh generated applies 
to both the thermal analysis and the creep analysis. 


INPUT! 



READ (5, 21) NPPE 
FORMAT (I 10) 


NPPE = number of pressure points per element 
This variable can be either 0, 1 or 3. If zero, the penalty 
function approach will automatically be used during the analysis. When 
1, the pressure will be constant throughout the element and compress- 
ibility will be satisfied only in an average sense in each element. 

When 3, the pressure is approximated by a linear function of x or y 
within each element. This last approach will insure satisfaction of 
the constraint of incompressibility everywhere within an element pro- 
vided the mid-side nodes are truly mid-side nodes of a triangular 
element (i.e., curved sides will prevent the constraint from being 
satisfied everywhere within the element— although NPPE = 3 can still 
be used for such elements and will more nearly satisfy this constraint 
than NPPE = 1) . 
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RPI = 0.5 
NDIVTH ° 11 
HDIVR o 4 



NDIVTH = 11 
HDIVR ■ if 



RPI ° 0.5 
NDIVTH = il 
NDIVR a 4 



Figure 1. Meshes generated by MESH1. 
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INPUT: 



READ (5, 3) RI, RO, RM, RPI 
FORMAT (4E10. 3) 


RI = inside radius 
RO ® outside radius 
RM = raid-radius 

RPI = multiple of it radians to be covered by the mesh. 

The use of those variable is illustrated in Figure 1. Note 
that RM is used to grade the element size in the radial direction. 
RPI * 2.0 will create a completely circular mesh. 


INPUT: 



1 PUNCH = 7 results in all mesh output data being read to 
TAPE7 for permanent storage 


INPUT: 



READ (5,1) NDIVTH, NDIVR 
FORMAT (2110) 


NDIVTH » number of mesh divisions in the theta direction 
NDIVR = number of mesh divisions in the radial direction 
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INPUT: 



READ (5, 10) NPBCI, XBCI, YBCI, TBCI 
READ (5, 11) CHI, CXI, CYI, TXI, TYI, TQI 
FORMAT (I 10, 3E10.3) 

FORMAT (6E10. 3) 


INPUT: 



READ (5, 10) NPBCO, XBCO, YBCO, TBCO 
READ (5, 11) CHO, CXO, CYO, TXO, TYO, TQO 
PORMATCHO, 3E10.3) 

FORMAT (6E10. 3) 

• 


These next four cards define boundary conditions along the inside 
radius and the outside radius. The vector components axe referred to 
the nodal point coordinates defined in Fig. 1. 

NPBC defines the type of boundary condition and the meaning 
of XBC, YBC, and TBC 

NPBC = ±1 both x' and y 1 components of velocity are unknown 

NPBC = ±2 x* component of velocity known to be XBC 

NPBC = ±3 y' component of velocity known to be YBC 

NPBC » ±4 x’ and y' components of velocity known to be 
XBC and YBC respectively 

NPBC < 0 temperature known to be TBC 

If XBC and YBC are not designated by NPBC to be known components 
of velocity, they can have one of two different meanings depending on 
the value of CX and CY. 

If CX = 0.0, then XBC equals known component of force in the x’ 
direction 

If CX 0.0, then XBC equals x' component of velocity external 
to boundary, to be used in viscous force boundary 
condition (see Eq. 10) 
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If CY « 0.0, then YBC equals known component of force in the y* 
direction 

If CY ^ 0.0, then YBC equals y» component of velocity external 
to boundary, to be used in viscous force boundary 
condition (see Eq. 10) 

CX “ coefficient of boundary viscosity in x* direction. It 
is the ratio of the boundary surface traction in the x' 
direction to the difference between XBC and the actual 
boundary velocity in the x* direction 

CY *» coefficient of boundary viscosity in y* direction. It 
is the ratio of the boundary surface traction in the y’ 
direction to the difference between YBC and the actual 
boundary velocity in the y' direction 

TX “ surface traction (farce/unit area) in the x' direction 

TY a surface traction in the y* direction 

When NPBC > 0, then TBC can have one of two different meanings 
depending on the value of CH. 

If CH * 0.0, then TBC is a point heat source 

If CH f 0,0, then TBC equals the ambient temperature used for a 
convective heat transfer boundary condition 
(see Eq. 24) 

CH = coefficient of convective heat transfer 
TQ = heat source per unit area 

C0SXXP = cosine between the local x' axis and the global x 
axis. The x' axis should always be defined so that 
this describes an angle in either the first or second 
quadrants, i.e,, the sine will always be assumed positive 


INPUT: 



RE AD (5,1) NUMBC 
FORMAT (110) 


NUMBC = number of nodal points which are to have individual 
boundary conditions specified. 
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INPUT*. 



READ (5,24) 11 1 NPBCC1I), COSXXPCin, XBC(I,1), 
YBC(Il ) , TBG (II) 

READ (5,27) CH(I1), CX(I1), CY(I1), TXCI1), 
TY(I1), TQ(II) 

PORMAT(2I10, F10.4, 3010.3) 

FORMAT C6E10. 3) 


As many pairs of the above r ;rds will bo read as specified by NUMBC. 

II * nodal point number 

NPBC(Il) = type of boundary conditions (see page 30) 

COSXXP(Il) = cosine between x and x' axes. See Fig. 1 
for example. The x' axis must always be 
defined so that the Sin(x ( x’) will be 
positive 

The remaining terms are defined on pages 31 and 31. 
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PR OGRAM MESH2--Input 

Program MESH2 is an automatic mesh generator. It is also through 
this program that the boundary conditions are specified. However, 
there are other opportunities to alter both the mesh and the boundary 
conditions. To understand the following input variables, one should 
be aware that two separate meshes are generated. The first mesh will 
always be contained within the second mesh and w*ll be used in CREEP, 
The second mesh will at least coincide with the first but can be 
extended beyond the first. It will be used in TEMP. 


INPUT: 



READ (5, 21J NPPE 
FORMAT (HO) 


NPPE = number of pressure points per element 
(See discussion under MESH1.) 


INPUT: 



READ (5,1) NUMLPS(l) , NUMLPS (2) 
FORMAT (21 10) 


NUMLPS (1) = number of loops in the finite element mesh for 
the plastic flow analysis 

NUMLPS (2) = number of loops in the finite element mesh for 
the thermal analysis, 
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In all cases, NUMLPS (2) >_ NUMLPS(l) , whether a thermal analysis is to 
be conducted or not, NUMLPS (1) should be specified as zero if only a 
thermal analysis is to be conducted. 


INPUT: 



READ (5, 3) XMIN, XMAX, YMIN, YMAX 
FORMAT (4E10. 3) 


These variables specify the minimum and maximum coordinates for 
plotting the mesh layout in an x-y Cartesian reference frame. These 
same variables are used for specifying coordinates for plotting 
results from program COUPLE. The exact use for these variables 
depends on the computer graphics available. 


INPUT: 



READ (5, 19) I PUNCH 
FORMAT(IIO) 


IPUNCH * 7 results in the nodal point coordinates and the 
nodal points associated with each element being 
read onto TAPE7. This output is not necessary for 
an analysis but is included in case it is desirable 
to have this data on punched cards or permanent file 
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THE FOLLOWING DATA WILL BE NEEDED FOB EACH LOOP IN THE ODDER PRESENTED *** 


INPUT: 



NDIV(I,1) * number of divisions on sides 1 and 3 of loop I 
NDIV(I,2) « number of divisions on sides 2 and 4 of loop I 
(See Appendix A for description of loops.) 

INPUT: 



This card indicates which sides of loop I are ’’-oined" to loops 
already formed. It is best explained by example. In the following 
figure, there are four groupings of numbers, each group representing 

/ 2 3 Z 1 0 0 4 2 

Side 1 Side 2 Side 3 Side 4 
the 4 sides of the current loop I. The card indicates that side 1 of 
loop I joins loop 2, side 3 (3rd side of loop 2) ; side 2 of loop I 
joins loop 3, side 1; side 3 of loop I does not join any loop already 
formed; and side 4 of loop I joins loop 4, side 2. The card only 
’’looks back” at loops already formed. It should not specify what 
loops not yet formed, i.e., loops > I, will connect to loop I, The 
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array has the following meaning: Loop I, side K, is joined to sido 

J0IN(I,J,2) of loop J0IN(I,J, 1) 

INPUT: 

/ 

READ (5, 11) (XC0R(I,J), YC0R(I,J), J ■ 1,8) 

FORMAT (2010. 3) 

i — — — — — w— i bum — — — — whim — ****** 

Each loop is a quadralatoral with parabolic sides. The geometry of the 
loop is defined by the eight points shown below. The eight data cards 
read in above define the x and y coordinates of these eight points, 
The first point read must always be a corner point. It, and the next 
two points, will define side 1 of the loop. The coordinates must be 
read counterclockwise, as shown: 



NUMBC = number of sides of loop I which have boundary 

conditions to be read. If a loop has zero sides 
which need boundary condition specified, NUMBC = 0 
will cause the next input to be skipped. 
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INPUT: 



f 

READ (5, 24) JI, LPBC(I,J), C0SLBC(I,J) 

XLBC(I,J), YLBC (I, Jl) , TLBC(I,J1) 

READ (5, 27) , CHBC(I,J), CXLBC(I,J1), CYLBC(I,J1), 
TXLBCfl, Jl) , TYLBC(I,J1), TQLBC(I,J1) 

FORMAT (110, 5E10.3) 

FORMAT (6E1 0.3) 



These two input cards are listed together since they are read in pairs 

for as many sides of loop I as specified by NUMBC. 

JI « side of loop for which the boundary 
conditions apply 

The remaining variables will be used to specify their corresponding 
boundary condition on all nodal points on side Jl. They correspond to 
NPBC, COSXXP, XBC, YBC, TBC, CM, CX, CY, TX, TY, and TQ which are 
defined on page 30. It should be noted that the corner nodal points 
will have the boundary conditions corresponding to the highest numbered 
side to which they are contiguous. These boundary points as well as 
other special points may be specified individually with the set of data 
cards which follow. 

THIS CONCLUDES THE DATA CARDS READ FOR EACH LOOP SPECIFIED BY NUMLPS(2 )*** 


INPUT: 



READ (5,1) NUMBC 
FORMAT (110) 


NUMBC = number of singular nodal points which have 
boundary conditions to be read. NUMBC = 0 
will cause the next input to be skipped. 
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INPUTS 


/ 


READ (5, 24), 11, NPBC (11), COSXXP (11), 
XBC (11) , YBC (11), TBC (11) 

READ (5, 27) , CU(U), CX(I1), GY(U), 

TX Cl 1) , TY(I1), TQ (XI) 

FORMAT (21 10, 4E10.3) 

FORMAT (6E1U. 3) 



As many of those pairs of cards will be road as specified by the 
previous value of NUMBC. II is the nodal point number which tho vari- 
ables, NPBC, COSXXP, XBC, YBC, TBC, CH, CX, CY, TX, TY, and TQ will be 
assigned. Explanation of those variables is given on page 30, 
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PROGRAM WAVE— Inpi 

Program WAVE takes the data from MESH and plans the tapo sections 
which will bo nocessary to store tho large finite clement matrices in 
programs CREEP and TEMP. 


INPUT! 

/ 


READ (5 >4) MSHCD (1) , MSHCD (2) 
FORMAT (21 10) 


MSHCD (1) ■ 1 identifies that program CREEP will bo called 
and hence tape segments should be planned 

MSHCD (2) ■ 1 identifies that program TEMP will be called 
and tapo segments should be planned for it 

Zero should bo specified when these programs are not to be 
called. 


The remaining READ statements in ;AVE will be read once for CREEP 
if MSHCD (1) = 1, and once for TEMP if MSHCD (2) a 1. 


INPUT: 



READ (5,4) MAXVOL, illMIN 
FORMAT (21 10) 


MAXVOL = the storage specified for the finite element stiff- 
ness matrices. For program CREEP it is the dimension 
of SKXX, SKXY, SiCYX and SKYY which must all be equal. 
For TEMP, it should bo one-half the sum of these 
dimensions. Actually, TEMP will use the storago of 
all four arrays for its stiffness matrix, however 
WAVE calculates the size of each tape segment based 
on half the storage available for this nonsymmetrical 
matrix 
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IBMIN « the minimum bandwidth to bo resorved In forming the 
K matrix. This normally should bo zero but if it 
happons that a mosh is so constructed that it will 
allow a very small bandwidth to bogin with and thon 
suddenly roquiro a largo bandwidth, this variable 
should bo sot equal to tho largor number. 


INPUT: 



READ (5, 4) 1READ 
FORMAT (110) 


Tho finite olemont matricos aro formed clement by element with each 
row, or equation, corresponding to a particular nodal point. In ordor 
to keep the bandwidth for these matrices as small as possiblo, it is 
necessary to choose the ordor in which the elements aro taken so that 
a minimum of rows, or equations, separate any two rows representing 
nodal points from the same element. The element numbering specified 
by MES1I1 provides the best ordor for this to occur. Therefore IREAD 
should be rood as zero which indicates the elements should be assembled 
in numerical order, and there is no need to read in a different ordor. 
However, when MESH2 is used, the element numbers may not indicate the 
best ordor to assemble the elements. In these cases, IREAD should be 
specified as 5. This designation will activate the next READ statement, 
otherwise it is skipped. 


INPUT: 



READ (5,11) (IORDER(IO,I) , 1=1, NUMEL) 
FORMAT (1015) 
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This is the ordor which tho elements will bo taken if other than 
thoir numerical order, See previous discussion on IREAD, Because 
two different moshos for CREEP and TEMP are used, IORDER must bo 
specified for each. A discussion on how to select a suitable order is 
treated in Appendix A. 
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PROGRAM COUPLE— Input 

Pry t ;cara COUPLE is the executive program which "couples" programs 
CREEP and TEMP, It also monitors the iterations between these programs 
to determine when steady- state conditions have been obtained as well 
as controlling the increments in a time dependent problem. All output 
is controlled by this same program. 


INPUT: 



READ (5, 12) NUMAT 
FORMAT (I 10) 


NUMAT = number of different materials to be identified 
in the analysis. Each element has a "material 
number." 


INPUT: 



READ (5, 22) (MAT (1) , 1 = 1, NUMEL) 
FORMAT (20 14) 


MAT (I) = material number for element I. If NUMAT = 1, 
this READ statement is by-passed and the array 
is initialized to 1. 


INPUT: 



READ (5 ,16) XNCPR, INCPU, INCPL 
FORMAT(3I10) 
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Those variablos stand for INCr oment to PRint, PUnch, and PLot 
rospectivoly. They represent the increments in the analysis that 
Subroutine PPP will print, punch and/or plot the output. The values 
read in at this time will be the first increment that these outputs will 
bo executed. Thereafter, they will bo executed at regular intervals 
(see next input). For a steady-state analysis, increments represent 
iterations between CREEP and TEMP, In a transient problem, increments 
represent steps in time. If it is not desirable to execute one of these 
routines before completion of the analysis, the corresponding INC-vaiue 
can be read in as zero (or less), or may be specified extremely large. 


INPUT: 



READ (MS) INTPR, INTPU, INTPL 
FORMAT (3210) 


These variables represent the interval, or number of increments, 
between each printing, punching, und plotting of the output. For 
example, if printing is called for at increment equal to INCPR, it will 
not be called for again until INCPR = INCPR + INTPR. This output 
is always called after the last increment of the analysis regardless 
of the values of INCPR, INCPU, and INCPL. If it is desired that one of 
these outputs not be used, even at that time, the corresponding 
INTerval number should be read in as a negative number. For instance, 
if it is not desired to have any punched card (or TAPE7) output, INTPU 
should be specified as -1. 
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INPUT: 



These five variables control the type of solution desired, be it 
a pure creep solution, pure thermal analysis, coupled steady-state, 
transient, etc. 

ITV = +1 will causo program CREEP to be called from COUPLE 

= -1 will cause program TEMP to be called. During a 

coupled problem, this variable oscillates back and 
forth between plus and minus 1. The value read in 
here will determine which of these two programs 
will be called first. In the case of a pure thermal 
analysis the value must be -1 whereas for a pure 
creep analysis the value must be +1. 

MOP stands for minus or plus. After either CREEP or 

TEMP is called by COUPLE, ITV is multiplied by MOP. 
Therefore, MOP = -1 will cause ITV to oscillate and 
produce a coupled analysis, whereas MOP = +1 will 
let ITV remain the value initially read in and 
hence will dictate either a pure thermal analysis 
or a pure creep analysis. 

TRANS is a variable specifying if a transient analysis is 
desired. TRANS =0.0 produces a steady-state 
analysis whereas TRANS = 1.0 produces a transient 
analysis. 

THETA is a variable that determines the type of numerical 
integration with respect to time that is used for 
the heat equation. 

THETA = 0.0 specifies a pure explicit integration 
THETA = 1.0 specifies a pure implicit integration 
THETA = 0.5 specifies a Crank-Nicolson integration 
It is recommended that THETA be taken between 0.5 
and 1.0. For a complete discussion of this variable 
see Crandall, 1956. 

INTEMP is a flag to indicate whether the temperatures in a 
coupled problem should be initialized at their 
steady-state pure-conduction values. If INTEMP = 1 
such an initialization will take place (it will 
involve a finite element solution in itself). If 
INTEMP = 0 no such initialization will take place. 
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LAGUEL stands for Lagrangian-Eulorian analysis. This type 
of analysis allows the mesh to follow the plastic 
flow in a transient analysis. It therefore corre- 
sponds to a control mass approach and is parti > j- 
larly convenient when the geometry of the problem is 
changing with the flow. LAGEUL * 1 implements such 
a solution procedure, LAGEUL = 0 will cause a pure 
Eulerian analysis to be made in which the elements 
are taken as divisions of space rather than divisions 
of material. Control volume analyses are made using 
this latter designations, 

IRZ is a variable which specifies whether the coordinate 
system is axisymmetric or not. IRZ = 1 indicates it 
is. IRZ = 0 will result in a plane flow analysis 
being made. In the case of an axisymmetric analysis, 
x corresponds to the r direction. The axis of 
symmetry is taken as the y(=z) axis. 


INPUT: 



READ (5, 26) TIMEM, MNI, DUMAX, DTMAX 
FORMAT (El 0.3, 110, 2E10.3) 


TIMEM = maximum time to be simulated during a transient 

analysis. For a steady-state analysis this variable 
can be set equal to zero. 

MNI = maximum number of increments. This variable limits 
the number of iterations during a coupled steady- 
state solution or the number of time steps in a 
transient analysis. 

DUMAX = the maximum displacement any nodal point is to be 
given during an increment of a LAGEUL analysis. 

DTMAX = the maximum value that is to be used for an increment 
of time during a transient analysis. The current 
version of the code allows only a constant value to 
be used for the increment of time. 
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INPUT: 

/~ 

READ (5,1) VECTL, CTEMP 
FORMAT (2E10. 3) 


Theso variables are used for plotting the output, and may be changed to 
fit the specific computer graphic equipment of the user. 

VECTL = the length of the arrow representing the maximum 

velocity of any nr dal point. It has the same units 
as are used for the nodal point coordinates. A 
good value to use is the shortest distance between 
ny two nodal points. 

CTEMP - the contour interval to be used for isotherm plots 
of the temperature field. All Isotherms will be 
multiples of this variable. 


INPUT: 



The code allows for the analysis of nonlinear creep behavior and 
nonlinear thermal behavior in addition to the nonlinearity resulting 
from the coupled problem. In these cases it is often desirable to 
iterate several times on either the creep equations or the thermal 
equations before returning to the other. ITMAXC specifies the maximum 
number of iterations to be made on the creep equations in such cases 
and ITMAXT represents the maximum number of iterations to be made on 
the temperature equations. 
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INPUT: 



Those six variables control how often an LDU decomposition of 


the finite element matrices should be made. 

INTLCU(l) * interval for LDU decompositions of creep matrix. 

This interval specifies the number of increments 
between CREEP and TEMP before a new decomposition 
is performed. 

INTLCU(2) = interval for LDU decompositions of creep matrix. 

This interval specifies the number of iterations 
within CREEP before a new decomposition is 
performed . 

INTLTU(l) = interval for LDU decompositions of temperature 
matrix. This interval specifies the number of 
increments between CREEP and TEMP before a new 
decomposition is performed. 

INTLTU(2) = interval of LDU decompositions of temperature 
matrix. This interval specifies the number of 
iterations within TEMP before a new decomposition 
is performed. 

LCU = 1 should be specified when the stiffness matrix 

for CREEP has not been saved from a previous run. 

= 0 should be specified when the LDU decomposition 
of the stiffness matrix in CREEP has been saved 
on a permanent file equated to TAPE2. 

LTU is the same as LCU but applies to the stiffness 

matrix in TEMP. If it has been saved on permanent 
file, it should be equated to TAPE4. 
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INPUT; 



READ (5, 31) DFCONV, DQCONV, DUCONV, DTCONV 
FORMAT (4E10. 3) 


These four numbers are convergence limits. They represent the 
minimum values for changes in 1) the right-hand side of the creep 
equation (DF = change in force) , 2) the right-hand side of the tempera- 
ture equation (DQ = change in heat source) , 3) the left-hand side of 
the creep equation (DU = change in velocity), and 4) the left-hand 
side of the temperature equation (DT = change in temperature) . In the 
current version of the program not all four limits are used. DFCONV 
specifies convergence for iterations within program CREEP such as 
needed for nonlinear flow laws. It is used in conjunction with ITMAXC 
to control the iterations within CREEP. DQCONV is currently not used, 
It would normally be used to determine convergence within program TEMP 
by comparing it to the total change in the calculated heat sources 
between each iteration. It is provided for future use whenever the 
need might arise; for the present, the variable may be read in as zero. 
DUCONV is al c n inoperative but allows the possibility of adding a con- 
vergence check by comparing the change in velocity between iterations 
of CREEP or between iterations of CREEP-TEMP when a steady-state analy- 
sis is being conducted. DTCONV is currently used to test convergence 
within TEMP and also convergence for steady-state analyses involving 
coupled behavior between CREEP and TEMP. This latter use is the most 
common use for this variable. It might be noted here that there are 
few needs for a nonlinear temperature analysis, therefore ITMAXT most 
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likely will bo specified as zero which will override any control 
DTCO.NV has within program TEMP* Therefore, DTCONV can be used 
exclusively for convergence tests on coupled thermal-mechanical analyses. 


INPUT: 



READ (5, 4) NPTS, NSEC 
FORMAT (21 10) 


The remaining INPUT cards allow final adjustment of either the 
mesh or values of the independent variables before execution of an 
analysis . 

NPTS » number of nodal point adjustments to be made 

NSEC * number of sections or groups of nodal points which 
are to be initialized to the same values 


INPUT: 



READ (5 ,35) II, XORD(Il), YORD(Il), UX(I1), UY(I1), UT(I1) 
FORMAT (IS, 2E10.3, 3E18.10) 


As many of these cards will be read as specified by NPTS. The FORMAT 
is the same as that used for TAPE7 output, hence it is at this point 
that output from a previous run may be read in as initial conditions 
for continued analysis, either transient or nonlinear. If punch cards 
are not to be used, it is best to write and read onto and off of 
TAPE7 directly without format. These cards and/or the next input 
can also be used to make changes in only a select number of nodal 




point values, for instance to insure convection will begin In a 
particular pattern. II equals the nodal point to which the values wil 
bo givon. 


INPUT: 



READ (5 ,9) JBGN, JEND, INCH, TEMPO, UXO, UYO 
FORMAT (31 10, 3E10.3) 


As many of these cards will be read as specified by NSEC. Each card 
will activate the following DO- loop. 

DO 185 J - JBGN, JEND, INCR 
UT(J) * TEMPO 
UX(J) = UXO 
UY (J) = UYO 


185 CONTINUE 
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Output Description 

All variables read into MANTLE are immediately "echoed" as output. 
Slosh marks (/) appearing within this echoed output indicate separation 
of input data cards. Each new lino of echoed data always represents 
a now input data card. Because theso variables have already been 
defined in the previous section they will not bo defined here. 


PROGRAM MESH1 or MESH2-- Output 



NUMNP a number of nodal points. This number is equal to 
NUMVP + NUMPP or NUMTP, whichever is greater 

NUMTP = number of temperature points, i.e., nodal points 
which specify temperatures in a thermal analysis 

NUMVP = number of velocity points, i,e, } nodal points which 
specify velocity in an analysis of creeping flow 

NUMPP = number of pressure points, i.e., nodal points which 
specify pressure in a creeping flow analysis 

NELMC ss number of elements specified for the creep analysis 

NELMT = number of elements specified for the temperature 
analysis 
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OUTPUT: 



r 

TEC 

rx 

TV 

TO 

CN 

cx 

cr 


i 900E+03 0. 

9, 

0. 

Ot 

0, 

0. 



> 900EI03 0. 

0. 

0, 

0, 

0. 

Oi 



1 900E+03 0. 

0. 

0. 

0. 

0. 

0. 



.900E+03 0, 

0. 

0. 

0. 

0. 

0. 



i 900E+03 0» 

0. 

0. 

0, 

0* 

0, 



• 900E+03 0. 

0, 

0, 

0. 


■ — ft_ " 



.9Q0E+03 0, 

0, 

0, 

0^_- 





i 900E+03 0, 

0, 

0. 






•900E+03 0. 

0. 







i 900E+03 0. 

0. 







i 900E+03 0, 
.900E+03 0. 
i 900EI03 0. ^ 
.900E+03 
*900F^^ 
i"X 








This is a complete listing of all nodal point data and parameters 
generated in program MESH. The only variables not already defined are 
the XORD and YORD arrays which are the x and y coordinates of the 


nodal point. . 
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OUTPUT: 



or 



This is a listing of the NP-array which gives the nodal point 
numbers associated with each element. The first six numbers are the 
nodal points associated with the velocities and/or temperatures and 
they are listed in a counterclockwise sense. The next one or three 
numbers appear only for elements used for a creep analysis for which 
NPPE / 0. They represent nodal points for pressure. 
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PROGRAM WAVE— Output 

Much of tho output from WAVE is not used by tho analyist and is 
providod more for tho benofit of those who wish to know how tho tape 
segments (blocks) are organized. This output is also helpful to deter- 
mine storage requirements to within critical valuos. Tho following 
output is printed once for tho creep analysis and once for tho thermal 
analysis if both are called for. 


OUTPUT: 


HBEQ If) LIOTX ICQHP IELEX HOVEX IEM PT 
IELEII) LIST! I) NPR< I) hOVE<I> INTO (If 


TAPE 

SEGMENT 

1 

KVOL 

EQUALS 

1345 


1 

35 

54 

24 

19 

0 

0 



I 

2 

3 

4 

5 

6 

7 

0 

9 

2 

42 

02 

122 

162 

202 

242 

202 

322 

124 

143 

203 

144 

204 

244 

1 

ai 

41 

363 

0 

A 

403 

443 

403 

401 

441 

5 

45 

05 

1/ 

TAPE 

SEGMENT 

2 

KVOL 

EQUALS 

1220 


2 

40 

50 

24 

10 

30 

26 



20 

21 

22 

23 

24 

25 

26 

27 

28 

243 

203 

204 

323 

363 

324 

364 

404 

403 

124 

145 

205 

164 

1 

01 

41 

161 

121 

345 

405 

445 

405 

7 

47 

07 

127 

167 

34 

35 

54 

36 

41 

42 

43 

44 

47 

33 

37 

30 

39 

40 

45 

46 

53 

54 

1 

2 

3 

4 

5 

9 

10 

11 

14 

31 

32 

33 

34 

35 

36 

37 

30 

39 

TAPE 

SEGMENT 

3 

KVOL 

EQUALS 

1260 


3 

40 

51 

23 

11 

26 

24 



30 

31 

32 

33 

34 

35 

36 

37 

1? 

\46 

244 

245 

205 

204 

325 

345 

126 j 


\ 

127 

120 

i 

01 

41 

141 

x 

' 


-^4? 

407 

447 

407 

7 








40 

41 , 

— 





10 

11 

12 

13 

14 

15 

16 

17 

362 

402 

442 

402 

3 

43 

4 

44 

141 

121 

241 

201 

243 

2B3 

323 

321 

125 

165 

205 

245 

205 

325 

204 



29 





443 

403 

444 

404 

5 45 

241 

206 

201 

321 

201 



40 <1? 50 

55 / 

15 1? / 

-10 


This output shows the organization of each tape segment. The 
first two lines of headings describe the variables which are listed 


for each segment. Each listing begins with the tape segment number 


and KVOL which is the actual volume needed to store the matrix for that 


particular segment. For efficient use of storage this value should be 
"approximately" the same for all segments with the exception of the 


last segment. 
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The first row of six variables aru defined by the first heading 
lino mentioned above. They are: 

NSEG « segment number which has already boon printed on the 
lino with KVOL, 

ID « tho maximum bandwidth for that particular segment, 

LISTX « the size of tho array LIST, This is tho number of 
equations in tnis tape segment. 

ICQMP = tho number of completed equations in tho current tape 
segment which are ready for tho LDU decomposition, 

IELEX a the size of tho array IELE. This represents the 
number of olomonts which will bo added during this 
segment in the formation of tho matrix, 

MOVEX - tho size of tho array MOVE. This is also tho size of 
the array INTO and represents the number of rows and 
columns of the old segment of the stiffness matrix 
which must be moved to be computable with the new 
segment. 

IEMPT = the number of empty rows of the current segment duo 

to the completed rows of the last segment having been 
read onto disc, 

The next several lines of output are the arrays indicated by the 
second heading line at the beginning of this output. 

IELE (I) = the elements added during this tape segment. The 
number will equal IELEX. 

LIST (I) = the nodal point numbers in this tape segment listed 
in the order they appear in the matrix. The 
number of entries equals LISTX. 

NPR(I) The printing of this array has been suppressed in 
the current version of the program. Its size is 
equal to the number of nodal points. Each entry 
has the following meaning: 

NPR(I) = 0 indicates nodal point I has not yet 

been brought into the stiffness matrix, 
i.e., the front has not reached it. 

NPR(I} t -1 indicates that nodal point I has been 
brought into the stiffness matrix and 
the decomposition has been completed 
and read onto low speed storage, i.e., 
tho front has passed it. 
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NPR(I) * somo positive integer value indicates that 
nodal point I is in tho current segment 
(or front) and is located in the LIST 
array at the value indicated, This corre- 
sponds to tho row of the matrix which 
corresponds to nodal point I. 

MOVE (I) and INTO (I) indicate that row and column 

MOVE (I) should be "moved" "into" 
row and column INTO (I). The 
size of each array equals MOVPX. 


OUTPUT: 



This output summarizes the previous output from WAVE and probably 
represents all that is of interest to the analyist. It is printed both 
for the creep mesh data and for the thermal mesh data. 

KMAX a maximum storage required for the stiffness matrix. 

For program CREEP this must be equal to or les* than 
the dimensions specified for SKXX, SKXY, SKYX, and 
SKYY. For program TEMP this must be evO&l to or 
less than one-half the sum of these four dimensions. 

IBMAX = the maximum bandwidth. 

NQMAX = the maximum number of equations that appui' 
any segment. 

NUMSEG = number of tape segments. 
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PROGRAM COUPLE, CREEP, TEMP—Output 

The remaining Output is the production data and the order in 
which it appears depends partly on the type of run being conducted. In 
the following discussion it is assumed that a coupled thermal-mechanical 
analysis has been made with program CREEP called first. 


OUTPUT; 


INCR 

0 

ITERC 

1 

DELU 

*101E-09 

DELF 

.121E4I7 

CP TIHE 
.274E402 

INCR 

0 

ITERT 

1 

DELT 

.22fiE+<M 

CP TIHE 
(4S1E+02 


INCR 

1 

ITERC 

1 

DELU 

.111E-0? 

DELF 
■ 133E417 

CP TIHE 
.650E402 

INCR 

1 

iteri; 

DELT 

.223E404 

CP TINE 
,824EK£_ 


INCR 

ITERC 

DELU 




This information is printed every increment of the analysis. 

INCR = the number of times control has been transferred back 
to program COUPLE from program CREEP and/or program 
TEMP. For a transient problem, this number will 
represent the number of increments of time taken 
since the beginning of the analysis. For a steady- 
state coupled analysis this number will represent the 
number of times both CREEP and TEMP have been called, 

ITERC - ITERT = the number of iterations performed in CREEP 
and TEMP respectively within a given INCR (without 
returning control to COUPLE) . 

DELU » the absolute maximum change in the value of any 

component of the velocity matrix as a result of the 
current iteration. 

DELF = the absolute value of the maximum change in the value 
found for any component of the force matrix as a 
result of the current iteration. 

DELT as the absolute value of the maximum change in temperature 
found at any nodal point between iterations in TEMP 
(printed only when there are iterations in TEMP) . 

CP TIME = the time on the central processor "clock." By 
subtracting from the current value, the previous 
value, one obtains the length of time taken for the 
current iteration. 
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OUTPUT: 



This is the general heading at the beginning o£ a complete listing 
of the production output variables. It will be printed as often as 
specified by INCPR (see Input discussion) . 

TIME EQUALS refers to current simulated time during a 

transient analysis. For steady-state analyses 
this will always be zero, 

NUMBER OF INCREMENTS EQUAL is self-explanatory. 

DELU and DELF summarize the final values of these two 
variables printed in the previous output. 

DELT = the absolute value of the maximum change found for 
any component of the temperature matrix during the 
last increment. 

DTIME a the magnitude of the increment of time between this 
current time and the previous increment. It has 
meaning only for a transient analysis. 

MAXIMUM DISPLACEMENT = maximum displacement given to any 
nodal point during the previous DTIME in a quasi 
Lagrangian analysis. 
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OUTPUT: 


H>f. 

XORb 

yord 

UX 

ur 

UT 

PTX 

FTY 

1 

- t 7i3/:. oj 

.347C+07 

•♦1 |4flfC"J4 

»234V2C-34 

.300000+04 

-.49313013 

♦IV444PI14 

2 

«943rt04 

»34:iH07 

-.2l604f>ll 

.34270012 

.30fl00r*Q4 

-.luitoma 

-.23432C414 

3 

t 107L+07 

♦ 310»:*07 

~.2Y/96f>»t 

.941113012 

.3UO0OUO4 

-.31477013 

.91430013 

4 

.none? 

♦ 307(;*07 

-»?r.v?w>u 

.13211011 

,30000004 

-,0rf?5E415 

.33334P+14 

a 

.2041:10/ 

♦2i)ino7 

-♦nitwic-n 

.96799013 

4300D0C4Q4 

«.2iAvintn 

-. 431171414 

A 

♦245T *07 

,?4*iti07 

♦ 3S7Mr-73 

-♦u/r/jr-n 

.30000004 

-.497114013 

.127431:414 

7 

,?nino7 

.204t*Q7 

+VAH04E-1? 

1312301) 

.300001.104 

-.?A792P+ir 

.169674:413 

0 

. JOVCI07 

♦ 13111*07 

.Mowr-n 

-.MAStiOU 

.300001.4 04 

-»7679H+l!i 

*51034013 

V 

.3IDE10/ 

♦ 10/1 *07 

♦vaj/nt-i? 

-♦av.-w-n 

.300001 404 

-.3933HHI9 

-ivn7iima 

10 

♦36 ir \ot 

.043006 

♦ 34Mt.£*l7 

-.21792 r-li 

.300001*404 

-.fii'.Aarur. 

-.54934013 

11 

.3471:10; 

-.70VJ>07 

♦ 3?rOJI£>34 

-♦640U/it.-14 

.30000C404 

,43M3+: + in 

-.lV3Alt:it4 

12 

.3411 *07 

-.MMIOA 

. 3413001? 

♦ JlfiMOU 

. 30000P404 

.imsirti* 

,354320)4 

13 

.sour »o/ 

-.10/4107 

«v6?2;r-i3 

,?V/7Qf>)l 

.300001:104 

.314011:410 

-.9140011 

14 

»30?t*07 

-♦lour»07 

* 1 J?04C'-13 

.sr.mc-u 

.auooot 104 

.U2724E4 IK 

-.313341+14 

IS 

♦ruino; 

-.204rto; 


.13100011 

.300002+04 

.21691015 

.431174:414 

16 

.24LE407 

-.24SC407 

-.33V10O13 

->339 tOC- 13 

.300001:404 

.6971+4013 

-.I2743P414 

17 

♦ro4Cio7 

-.2011107 

-♦13J24E-H 

-.960 02012 

.30000P+04 

.24792015 


Itt 

.t5lIC*0? 

-»309r407 

-.;*a634F'U 

Mllaottoll 

,30000004 

.76770C+I5 


IV 

♦ 107007 

-.330E+0/ 

-.292V2E-11 

«.9&mOl2 

.aoooono4 

. 33531+1: 


20 

♦ s43[rn'6 

-.343E+07 

-*2i7rir-ii 

-.34314012 

.30000C404 



21 

-.712E-QJ 

-.347E40? 

-.M3D7E-M 

.13001 P-23 

•30000E+04 



22 

-.B41EI04 

-♦343CI07 

.31035E-U 

-.34139012 

.3C0OCE4O4 



23 

1Q7C+07 

-.33*^407 

.29770011 

-.96729012 

.30000004 



24 

-.ISBE+O? 

-♦309C407 

.269)4011 

-.132040)1 

.30000E + P>^ 



33 

-♦204E+Q7 

-.2DIE+07 

» J 31 BQOl 1 

-.9Q74CO 12 

»30000J>^ 



26 

-i24*r+07 

-.243CIQ7 

-.33910013 

.33910013 




27 

-.MUTtO? 

-♦204E407 

-.96003012 

.133241-11 




2B 

-.307F4Q7 

-.irCF407 

-* 1306101 1 

.25434011 . 




SJi* 

~i330r.*07 

107E4O7 

-.9&I77C-12 

.27272E-L>^ 




x 

-.3I3E+07 

-.3431:406 

-.34514012 






SO4/C107 

.142E-02 

-.31232034 






^-<*07 

.043C40A 

-.341300*1^- 






This is tho listing of all nodal point production variables. 


Variables not yet defined are: 

UX = the x-component of the nodal point velocity in global 
coordinates . 

UY = the y-component of the nodal point velocity in global 
coordinates , 

UT = nodal point temperature. 


FTX, FTY are the x' and y' components of the reaction 
force acting on the nodal point. They are given 
in terms of the local coordinate axis defined 
by the direction cosine COSXXP(I). They are 
obtained by solving: 


where 


{FT} = 
F 


[K]{U} - {F g } - {F e > - {F a > 


S 


! the gravitational force 

F = the elastic force, and 

e 

F = the applied forces. 

3. 

Hence the FT forces arise only from reaction forces and should be zero 


except at constrained points. 


Because these values depend on the 


velocity vector, at least two iterations or increments must have been 
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made for thorn to have any meaning. For a nonlinear problem, those 
forces are in error until convergence has been reached. For this reason 
they sorve as a convenient measure of convergence and indicate where 
the error is greatest. 


OUTPUT: 



These remaining nodal points (which have the same numbering as 
some nodal points for temperature when the temperature mesh is larger 
than the creep mesh) are used for specifying pressure. For ease in 
programming,, each node has two pressures associated with it in the 
same manner as the components of velocity are specified. Both PRS 
and QRS are pressures associated with two different elements, designated 
by XPQ. Hence, in the example above, nodal point 30 is used to specify 
the pressure in elements 5 and 6. PRS is the pressure in element 5 
and QRS is the pressure in element 6. For average incompressibility 
(NPPE = 1) there is only one pressure specified per element and there- 
fore it is assumed constant throughout. For complete incompressibility 
(NPPE = 3), there are three pressures per element. Their values should 
be considered as located at the midside nodal points of the side 
connecting the elements designated by IPQ. 




61 


OUTPUT; 


f-H 

muii 

niuxxi 

ClllYYI 

nmxri 

nwxxs 

0IOYY2 

DMIXt? 

ninXKS 

stam 

810X73 

1 

tov 

, .3/ OF Mlfl 

■**421 noil 

iW/no; 

-t 104C 1 07 

* 1061: lOV 

-»V;’?E107 

• »ir?nov 

*n?noj 

tdJUOfl 

5 

• ansionon 

♦ 434O0N 

-»4S/non 

-.141000 

->4iur ion 

*4;/r ron 

.irnnon 

,7nnoa 

*,67/>nOR 

-»lfT6ri07 

3 

» J'iJAW 10(1 

*4?vuoti 

-, 4 imi ion 

,J/M“K >6 

i.iur ion 

-*4.tu< ion 

,2601 10/ 

",476U0« 

,4Vflt‘*0f> 

*3l»rM10 

4 

tWitM mu 

,Htvr MJfi 

km 

*4/01 1011 

-.rwut ion 

,2v /r Mm 

,2i7rion 

-.f.nuioo 

,614110(1 

,2/.inon 

3 

♦ JlHVOtlOU 

iL’OVUOU 

-.17IU0U 

,401,0 08 

4 lOiC 107 

-,i)4do? 

,2176107 

-.204E108 

.717CMIB 

il94Clim 

6 

*S0344CICU 

N.nunoy 

♦ u.t,ckxi 

,177000 

.lOEfOn 

-,220000 

,?cjoon 

-,M3EI0? 

,?12E407 

oiotton 

7 

,40017000 

-,3ivoorj 

,304000 

.WC1CD 

-.377000 

436*000 

, 22 anon 

4 1ICCI0U 

-4 170000 

", 40ftl07 

0 

«9U6firtOD 

-.rmrica 

•SflACICB 

iVPVEIO? 

* 04 ^ic ion 

-.nwrion 

»,4/4( f 0(1 

.6270011 

- » A4 1 E 1 00 

-»610E 100 

9 


"iSft-IOWI 

♦ 36M? 101) 

« •f J ; , oc ion 

,?/;'Udo 

-,YV4t)C(f 

-.nut 107 

.u^r.t i no 

-1077(100 

-noitiov 

10 

umsocv 

-»*7cnon 

.Drsooa 

-»43 r .CIQfl 

.avooon 

-,0S9dQH 

-.nwnoa 

-.JiKfOfl 

.377E100 

-,731E100 

n 

,1237700? 

.377E1O0 

-4340000 

-.426000 

.447O0U 

-»416CM>0 

-.j/roou 

,0461100 

-,703C100 

-•CMT*rS£. 

12 

,11937007 

*337000 

-.C17O0& 

-,172000 

.711000 

-.Anvooo 

471007 

,0070011 

"♦760C10«^ 


U 

,13333007 

,634001) 

-,43QElQa 

-*710007 

,?niooa 

-.777000 

.337000 

,nooc(OQ 


14 

.4**94008 

.200000 

-.sooEtoa 

.7400 07 

*304000 

-,304000 

,274E108 

♦ 413EI0I1 



13 

,5Vo.f400n 

,27*000 


tiotctoa 

.224000 

-.23OE1O0 

*170000 

4 307^^ 



n 

. I3JI1&CI0V 

,331000 

-.33*000 

*40:000 

,71)7E»07 

-liOOCICO 

,701000 




i? 

472337000 

,13/0011 

-» (740011 

,4270 00 

-.2000 00 

* 174F10D 

,310 Ut>y 




m 

.726370011 

<2ivci07 

-♦300E1O7 

,324000 

-.3Q4tfCD 

, root too 

»440>"^ 




IV 

,a3??DnaG 

.402007 

-.474007 

♦310E4O0 

» nor. loo 

-.120000 





so 

,63273008 

-.161000 

* 147000 

->2361*10? 

-.471000 

.MtOOP. 





21 

♦ 531*0000 

.33OE1O0 

-♦2V1O08 

»663E(0ff 

*4?4E40n 

-,443F>^ 





>03 

,7ia-1*,EfO£l 

»S2 toon 

-,4u?r»on 

* 37XC* 041 

♦ 377 000 






\ 

,10447007 

.40001) 

-,410000 

4270000 

.4770 0*, 







Simictw 

»34?ooo 

-.349O0B 

-.isvrm 








^iyocu 

•rvucion 

-,304EI0tt 

-.lOEIjM 









.SiflEIOQ 

-,27dEt0Q 









This figure shows the production output for each element. 


SIGH = the effective stress as calculated at the center 
quadrature point 

SIGXX, SIGYY, and SIGXY 1,2, and 3 are the three components 
of the stress tensor at their separate quadrature 
points. (See Appendix for designation of all 
quadrature points.) 
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Error Messages 

Most of the orior messages in MANTLE are self-explanatory. Those 
which are not solf-oxpinnatory or which require special action to 
correct, will now bo discussed. 



This message results from there not being sufficient room to 
rearrange the stiffness matrix from its orientation specified by the 
LIST array of the previous segment to the orientation needed for the 
segment under consideration. This causes an automatic STOP of MANTLE. 
The two LIST arrays printed are difficult to interpret and will not 
be discussed here. However* three pursuits to alleviate this problem 
can be made. 

1, Increase KVOL (and correspondingly the dimension of 
C0MM0N/C7/). 

When this is possible it is usually the easiest and best approach. 
If this is not possible, decreasing KVOL can sometimes remedy the 
situation. This occurs simply because any change in KVOL results 
in an entirely new arrangement of the. tape segments. It often 
happens that one value of KVOL will work when a KVOL either higher 
or lower will not. KVOL can be made as low as zero and still 
produce legitimate tape segments. 
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2, Chock to soo if tho first tapo segment has a bandwidth (*!;,) 
significantly lower than tho bandwidth of tho last segment 
printed. 

If so, increase IBMIN to equal tho largest IB so far printed. 

3. Specify a new order in which the elements are taken, to 
produce a smaller "front" (see Appendix for more explanation) . 

If all three methods have been exhausted the finite element matrix is 
too large for the computer. Tho only alternative, in this case, is 
to reduce the size of tho matrix. This can bo done by reducing the 
number of elements, or for a creep analysis, reducing NPPE to zero 
(penalty method) . 


OUTPUT: 



This is a warning message and will not cause the run to stop. It 
occurs when KVOL is greater than the specified MAXVOL. It is important 
to note that this is not a check on the actual dimensions used in 
COMMON/C7/. It can, in fact, be desirable to specify MAXVOL less than 
the available storage in order that other arrays associated with the 
size of tape segments remain within their limits. If KMAX does exceed 
the actual dimensions of C0MM0N/C7/, simply reduce the specified MAXVOL. 
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This message indicates that in lifting the order of the elements 
as data, the user has made an INPUT error. Automatic stop. 



This occurs when the determinate of the Jacobian matrix is found 
to be negative during the Gaussian quadrature. This will be the case 
anytime the nodal point numbers in the NP array are not listed in a 
counterclockwise order as well as when an element has a reentrant side 
of significant magnitude. These conditions are most likely to occur 
in a quasi-Lagrangian analysis when the elements have become extremely 
distorted due to the flow. Automatic stop. 
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This occurs during the Gauss elimination whenever a diagonal term 
is found to be zoro at the time that zeros are being placed beneath it. 
Strictly speaking this should novor occur and most likely means a 
major error in the compiled program. It can possibly arise if a given 
element is overly constrained, i.e., most nodal points given known 
velocities. Automatic stop. 



This occurs anytime that 200 or more contours are to be plotted 
in one element. This would normally represent an error in the analysis 
or a poor choice of contour interval. In either case the computer 
time to perform the plotting is often large and therefore an automatic 
stop is called for. 
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User Subroutines and Functions 

There 1 b always a balance between selecting a very general program 
where a wide variety of problems may be solved by entering sufficient 
Information as data, and selecting a very dedicated program designed 
for a select group of problems. The advantages of the general program 
are rather obvious but it haB the disadvantage of requiring a large 
amount of input data and It is usually much larger than a dedicated 
program. On the other hand, the dedicated program has the obvious 
disadvantage of being limited in its scope but has the advantage of 
being operable with a minimum amount of input, usually shorter run 
times, and less total storage. One method which incorporates the 
advantages of both methods is to make available to users the ability to 
write their own dedicated subroutines. This approach was used in the 
design of MANTLE und these subroutines will be discussed here. The 
user will be able to bring most of the pertinent information into these 
subprograms through the common blocks. That information not contained 
in common will bo brought in through arguments. 

FUNCTION GAMX (TEMPK, XK, YK, MJ) 

FUNCTION GAMY (TEMPK, XK, YK, MJ) 

These two functions represent body forces (a force per unit 
volume) in the x and y directions. They are called by CREEP during 
the matrix assembly at each quadrature point so that Gaussain quadrature 
can be used to integrate the total body force. This takes place in a 
J-DO-loop over each element and a K-DO-loop over the quadrature points. 
TEMPK = temperature at quadrature point K 
XK, YK e x and y coordinates of quadrature point K 
MJ = material number of element J 
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FUNCTION G (TEMPK, MJ) 

G 1b the shear modulus of material MJ at temperature TEMPK, This 
value Is needed for a steady-state elasto-vlsco-plastic analysis, 

There Is much research still to be done with regard to these analyses 
and the user should proceed with caution. To by-pass this analysis In 
CREEP (i.e., produce a purely viscous analysis) G should be set equal 
to a negative number, e.g,, 

G « -1.0 

This has the same effect as using a very large G but with much less 
computer cost. The use of G Is limited to steady-state analyses. 


SUBROUTINE VISC (VS, VT. PENLTY, NPPE, EPSII, TEMPK. XK, YK, MJ) 
Rather than a FUNCTION statement for viscosity, a subroutine was 
chosen so that more than one value could be returned to CREEP. These 
values are VS, VT and PENLTY and are defined below: 

VS * the viscosity based on a secant modulus. That is, VS is the 
value of Vi at quadrature point K to be used in 


and is, in general, a function of the material, MJ; the effective strain 
rate, EPSII; the coordinates XK, YK, and the temperature, TEMPK. When 
p is a function of the strain rate it is possible to use a viscosity 
based on the tangent of the stress-strain rate curve. That is, the 
tangent viscosity, VT, defined as 
d Ty - 2P T d Gij . 


The third variable, PENLTY, 1 b the penalty function. Its value is zero 
unless the number of pressure points per element, NPPE, is zero. In 
this case it is set equal to 1000.0 * VS, i.e., a large number in 
comparison to the viscosity. PENLTY appears as X in the equations 
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°ij ’ 2ue ij * Xe kk s ij 
da ij * 2 *'r‘ lE ij * *t de kk )s ij ’ 

It is scon therefore, that by oquating p and X to the elastic 
constants G and X it is possible to use CREEP for the analysis of 
olastic solids. Because six nodal point elements are used, full inte- 
gration is used for the penalty function rathor than a reduced 
integration. 

Two words of caution are necessary. In a viscous flow analysis, 
a small X doos not roprosent a compressible flow analysis. Wien 
a small X is used in an olastic analysis, the stresses calculated in 
MANTLE roprosent neither the total stress nor the stress deviator. 

FUNCTION RHOQU, 1ELEJ) 

FUNCTION CPCMJ, J ~ LEJ) 

FUNCT ION RKXCMJ, IELEJ) 

FUNC'i ION RKYCMJ, IELEJ) 


These four functions represent the density, specific heat, and 
thermal conductivity in the x and y direction respectively for 
material MJ, These parameters are assumed constant within each element 
and hence are called only once for each element (in contrast to VISC 
for example, which is called at each quadrature point). The element 
number is also made an argument in case it is desirable to know location, 
temperature or other information which can be obtained by knowing the 
element rather than the material. In these cases it is only necessary 
to include the appropriate COMMON blocks in the function programs. 
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SUBROUTINE MSHADJ 

This subroutine is called near tho end of both MESH programs and 
is provided in order that tins analyst can make adjustments in the 
finite element mesh, Those adjustments can tako a number of forms, 
tho most common being slight adjustments of boundary coordinates to 
fit unusual shapes. Other adjustments could bo specification of 
boundary conditions along curved boundaries which vary according to 
some specified function, changing a direction of one or more diagonals 
in the mesh, and adding additional elements. 

SUBROUTINE BNDRV 

SUBROUTINE BNDRYC 

SUBROUTINE BNDRYT 

These three routines are called during COUPLE, CREEP and TEMP 
respectively. BNDRY is called only once and can be used to make any 
adjustments necessary for a specific run. The next two are called 
every increment and are intended for the purpose of altering boundary 
conditions during iterative analysis, but can be used for other 
purposes as well. As an example, boundary conditions could be altered 
due to material coming into contact with other boundaries during the 
flow analysis. Free surface adjustments for steady-state analysis also 
are made within these subroutines. BNDRYC is called at the end of each 
iteration within CREEP, BNDRYT is called at the beginning of each 


iteration in TEMP. 
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SUBROUTINE STIFF (IEhHJ, ITV) 

This subroutine is called immediately after the stiffness matrix 
has been formed and before it has been assembled into the global stiff 
ness matrix. The element number is an argument as well as ITV which 
indicates whether a creep analysis or a temperature analysis is in 
progress. 

The purpose of this subroutine is to provide the user the 
opportunity to make unusual changes in the matrix before assembly. As 
examples, a penalty function could be used to force two or more compo- 
nents of velocity to be the same. Also, alterations in the matrix 
could be made to simulate a failure plane across the element. The use 
of this routine is limited to those users who are fairly familiar with 
the finite element method of approximation. 
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ILLUSTRATIVE EXAMPLES 

In this section soveral problems are presented to help the user 
establish some valid test runs. The examples are limited to mantle 
convection, however the user may consult the reference to obtain 
additional examples. 

Example 1, Pure Conduction 

This problem can be used to check SUBROUTINE TEMP. Figure 2 
illustrates the problem which is one of steady-state, pure conduction 
in an infinite cylinder. The dimensions used coincide with the earth's 
mantle. Bet... the inside and outside temperatures are constant and 
there is no convection. Table I compares the exact solution with that 
obtained by MANTLE. Any row of nodal points in the radial direction, 
such as emphasized in Figure 2, can be used for comparison. 


Inside: 


r = 0.347 

x 10 7 m 

T = 

900.0 °C 

Outside: 


r = 0.637 

x 10 7 m 

T = 

3^0.0 “C 



Figure 2 
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TABLE X. Solution to Example Problem I 




TEMPERATURE 

Nodal 

Point 

Radial 
Distance 
x 10“ 7 

MANTLE 

Solution 

Exact 

Solution 

1 

0.347 

900.00 

900.00 

41 

0.371 

833.54 

833.94 

81 

0.395 

771.26 

772,02 

121 

0.419 

712.67 

713.76 

161 

0.444 

657.36 

656.52 

201 

0.468 

604.97 

604.52 

241 

0.492 

555.22 

555.12 

281 

0,516 

507.85 

508.08 

321 

0.540 

462,64 

463.18 

361 

0.565 

419.41 

418.47 

401 

0.589 

377.98 

377.38 

441 

0.613 

338.23 

337.93 

481 

0.637 

300.00 

300.00 


Exact Solution: 

Temperature = -977.744 r + 15775.093 
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E xample 2. Thick-Walled Cylinder 

This problem can be used to check SUBROUTINE CREEP. Figure 3 
illustrates the problem which is a thick-walled cylinder under an 
external traction. Again, the dimensions are those of the earth's 
mantle and the element layout is the same at that used in the previous 
example. Table II compares the exact solution with that obtained by 
MANTLE, where the velocity is in the radial direction. 



Viscosity - 1.0 x 10 N*S/m 

7 

Inside radius = 0.347 x 10 m 

7 

Outside radius = 0.637 x 10 m 


Surfaces Traction 
T x = 700 N/m 2 


Figure J> 



74 


TABLE II. Solution to Example Problem 2 


RADIAL VELOCITY 


Nodal 

Point 

Radial 
Distance 
X 10- 7 

MANTLE 
Solution 
x 10" 4 

Exact 
Solution 
x 10“ 4 

1 

0.347 

0.17255 

0.17270 

41 

0.371 

0.16139 

0.16152 

81 

0.395 

0.15157 

0.15171 

121 

0.419 

0.14288 

0.14302 

161 

0.444 

0.13513 

0.13497 

201 

0.468 

0.12819 

0.12805 

241 

0.492 

0.12193 

0. 121t! J 

281 

0.516 

0.11625 

0.11613 

321 

0.540 

0.11108 

0.11097 

361 

0.565 

0.10636 

0.10603 

401 

0 . E&9 

0.10202 

0.10174 

441 

O f ft 1 — 

0.09803 

0.09776 

481 

0.637 

0.09433 

0.09407 


Exact Solution: 
velocity = 


General Solution: 


Velocity = 


'r.r (T 

X 0^ 0 


T.) 

i' 


r. 

l 


l 

r 


where 


r^ = inside radius 

r = outside radius 
o 

T. = inside surface traction 
x 

T « outside surface traction 
o 
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Example 3, Simple Convection 

This example as well as those that follow, deal with convection 
within the mantle. With the exception of the final example, it has been 
assumed that convection extends the entire depth of the mantle. Al- 
though there is some evidence that this might be the case, there is 
no Attempt in this report to try to prove or disprove that theory. The 
choice was made simply to provide some interesting examples with a 
minimum of computer cost. If the upper mantle is modeled, and an 
entire cross section of the earth is taken, many more elements are 
needed to provide a reasonable element aspect ratio than are needed 
for modeling the full mantle. 

The parameters used for these analyses (see Table ^.11) represent 
typical values as reported in the literature, the exception to this 
being the viscosity. In order for the iterative process used in the 
analysis to converge, the number of elements must increase with the 
Rayleigh number. For this reason, the Rayleigh number was kept near 
its critical value, which for full mantle convection requires a viscos- 
ity larger than that normally assumed for the mantle, However, the 
examples do indicate the versatility of the code, and are of interest 
in themselves. 

This example assumes a constant viscosity within the mantle equal 

25 2 

to 10 N*S/m . The fluid is allowed to slip at both the inside and 
outside boundaries. In order to decrease computer cost, a slightly 
coarser mesh was used than in the previous examples and is shown in 
Figure 4. In all examples, NPPE = 0, which results in the penalty 
function approach. There are, therefore, a total of 440 nodal points, 
giving 880 degrees of freedom for the velocity field. 
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C 


C 


c 


c 


c 


c 


c 


c 


FUNCTION GAMX<TEMPK»XKr YKrMJ) 

F«-0*74*TEMPK 

R=SGRT<XK*#2+YK*#2> 

GAMX«-F#XK/R 

RETURN 

END 


Dimension? 
pga = [kg/m 2 .g 2 °C] 
R, XK, YK [m] 
TEMPK [°C] 


FUNCTION GAMY < TEMPK » XK > YK r M J ) 

F 33 -0*74*TEMPK 

RESORT < XK*#2+YK**2 ) 

GAMY»-F*YK/R 

RETURN 

END 

FUNCTION G < TEMPK » MI) 

G=~l *0 
RETURN 
END 


GAMX, GAMY [kg/m 2 *s 2 
* N/m 3 ] 


G [N/m 2 ] 


FUNCTION RHO(MJrlELEJ) 

RHQ 33 ♦ 370E+04 RHO [kg/m 3 ] 

RETURN 

END 


FUNCTION CP(MJrlELEJ) 

CP=1.2E+03 CP [J/kg*°C] 

RETURN 

END 


FUNCTION RKX<MJrIELEJ> 

RKX*6.66 

RETURN 

END 


RKX, RKY [J/m* s* °CJ 


FUNCTION RKY<MJ»IELEJ> 

RKY® 3 <5 ♦ 66 

RETURN 

END 


SUBROUTINE UISC ( VS»UT» PENLTY r NPPE *EFSII* TEMPK r XK r YK fMJ) 
VS»5*0E+24 

UT-5 *0E+24 _ 

PENLTY«0*0 VS, VT [N*s/m ] 

IF<NPPE*EQ»0> PENLTY 33 1000. 0#VS 

RETURN 

END 


TABLE III. Parameters for Example Problems 3-7 
Excluding Viscosity 



77 


The stoady state isotherms are shown in Figure 5, Iterations 
were begun using a constant temperature field with four temperature 
"spikes" evenly spaced. After five iterations, the largest change in 
any nodal point temperature from the previous iterations was ld.5°C. 
The largest surface velocity was found to be 1.793 x 10"** m/s or 
0.05655 cm/ year. If the Rayleigh number is defined as 



where d = depth of mantle, i.hen, for this example R = 3610 or about 
5.5 the critical value. 




7 

Inside radius = 0.347 x 10 m 

7 

Outside radius = 0.637 x 10 m 


Figure 4. Mesh for Ext..nples 3-6 
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Example 4. Low Viscosity Zo ne 

This example was chosen to illustrate the use of a variable 
viscosity. Although viscosity can be specified as a function of temper- 
ature, position, and/or strain- rate, only position (in this ;ase, depth) 

was used for this examplo. The viscosity of the mantle at a depth 

25 2 

greator than 700 km was taken as 10 N*S/m as in the previous examplo. 
In the upper mantle, depths loss than 700 km, the viscosity was taken 
as 10 23 N*S/m 2 . 

The problem was initialized with the temperature field from the 
previous examplo. As the iterations continued, the eight-cell pattern 
broke into a four-cell pattern. For this particular case, steady-state 
was not reached, Figure 6 illustrates the results found on the 13th 
iteration. During the next iteration, the temperature became unstable 
and the solutions begin to diverge. However, it is believed that the 
overall flow pattern is reasonably correct. It is interesting to note 
that reducing the viscosity in the upper mantle resulted in larger con- 
vection cells. It is not unlikely that a finer mesh in the upper mantle 
would have resulted in smaller convection cells in that region, giving 
rise to a two-scale convection pattern. 
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Example 5. Axisymmctric Convoction 

Tho previous two oxamplos havo assumed plane flow. This oxample 

shows what effect tho constraint of piano flow can havo on the 

analysis. Only half of tho mantle cross section, as shown in Figure 7, 

25 2 

was usod, Tho viscosity was assumed to bo a constant 10 N*S/m . Tho 
problem was first analyzed as one of plane flow. The isotherms for 
this are shown in Figure 7a, Next, the flow was analyzed as axisym- 
motric. Those results are shown in Figure 7b. Again, the four-cell 
(eight-coil full cross section) convection pattern broke into a two- 
cell (four-cell full cross section) pattern. Both results were obtained 
to a convergence limit of less than 50°C between iterations. The 
maximum surface velocity found for the two cases was essentially the 
same at 0,06 cm/year. 
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Example 6. Crustal Platos 

This example was conducted to illustrate one use for SUBROUTINE 
STIFF, as tfell as to examine possible effects on mantle convection due 
to rigid crustal plates. The platos wore simulated by requiring that 
groups of nodal points on the surface have the same tangential velocity. 
Once again, plane flow was assumed, hence the "plates" are on the out- 

25 2 

side of an infinite cylinder. The viscosity was taken as 0,5 x 10 N»S/m . 
The steady-staje solution with no plates (similar to Example Problem 3) 
is shown in Figure 8, Figure & shows the plates that were first assumed 
and the isotherms and streamlines resulting from the analysis. Clearly, 
the additional constraint imposed by the plates results in an entirely 
different flow pattern. For this four-plate example, the flow was 
counterclockwise at all surface points. In an effort to eliminate 
this characteristic and thus create a more interesting flow pattern, 
one plate was removed. Figure 10 illustrates the new plate configura- 
tion with the resulting isotherms and streamlines. 

For both the four-plate analysis and the three -plate analysis, 
t, ,'Vergence was not obtained. Although the iterative process did not 
diverge, large changes in temperature (and thus flow patterns) resulted 
betweer each iteration. The figures shown are simply typical of thosj 
that were found. It is not fully under tood why convergence was not 
obtained but it is assumed that additional elements would help to 
alleviate the oscillations. One could study these oscillations by 
terminating the iterative procedure and beginning a transient analysis. 

This should eventually lead to a steady state condition. 



Contour interval = 300°C 


Figure 8 


Steady State Isotherms for 
Example 6, No Plates 
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Plate velocities in m/s 

Figure 10a. Streamline fox Three-Plate Analysis-- 
Example Problem 6 
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The surface constraints were imposed using the penalty function 
method. The tangential components of two points n and m are the 
same if 

UY (m) - UY(n) = 0 . 

This constraint can be incorporated into the variational principles by 
adding to the original functional term 

A (UY (m) - UY(n)) 2 

where A is assigned a large value. Perhaps an easier to understand 
explanation of the method is to consider an additional one-dimensional 
element which connects 'the two nodal points being considered. Figure 11 
illustrates such an element as a connecting link between the two points 


A 



Figure 11 
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with a viscosity A. Tho stiffness for this one-dimensional element 
is expressed by 

'X -X-j UY (m) FY (m) 

K 

.-X Xj UY(n) FY (n) . 

If X becomes large and there are no corresponding large terms added 
to the right-hand side, it is clear that the desired constraint will 
be realized. The small 2 ;i 2 stiffness matrix is simply added to the 
element stiffness in SUBROUTINE STIFF. 

One other aspect of the code is illustrated in this example. 

The nodal point boundary conditions specified in the array NPBC (I) 
will be unchanged by even additions of ±10, provided the sign is not 
changed. That is, the following specifications will accomplish the 
same nodal point boundary conditions: 

NPBC = 4 and NPBC = 14 and NPBC = 24 

or 

NPBC = -2 and NPBC = -12 and NPBC = -32 

or 

NPBC a 3 and NPBC = 43 and NPBC = 123, etc. 

This allows the user to specify special boundary conditions to be used 
in any of the subroutines provided. In the present case, any two 
adjacent nodal points in the same element with an absolute value of 
NPBC greater than 10 were taken to have the same y ' -component of 
velocity. 
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SUBROUTINE STIFF < IELEJfITV) 

COMMON/C1/ 

1 XORD< 550) rYORD< 550>rXBC< 550>.YBC< 550>rTBC< 550>r 

2 CX< 550) rCY ( 550 > r CH ( 550)rTX< 550>rTY< S50>rT0< 550>r 

3 COSXXP< 550) »NPBC< 550>rNP< 300r 6> 

COMMON 

1 TXX<9r9) rTXY(9r9) rTYX<9r9 ) fTYY<9i 9> r 

1 SXX<9r 9) »SXY<9r9) r SYX<9r9> r SYY<9r 9 ) , 

1 SPX<6r3> rSPY<6r3) r 

1 SIGXXJ<3> rSIGYYJ<3> rSIGXYJ<3) rSIGTHJ<3) , 

1 DNCJDX < 3 ) r BNGDY < 3 ) r 

1 RJAC<2r2)rRJACI<2r2>rDNDX<6> rDNDY<4> 

IF < ITV«EQ* -I ) RETURN 
BO 500 l~lr6 
NPI«NP(IELEJfI) 

NBCI«IABS<NPBC<NPI>> 

IF<NBCI «LT* 10) GO TO 500 
BIG®SYY ( I f I ) *10000*0 
DO 400 J«If4 
IF(I.EQ.J) GO TO 400 
NPJ=NP(IELEJfJ) 

NBCJ=»IABS(NPBC<NPJ) ) 

IF<NBCJ*L.T ♦ 10) GO TO 400 
TYY<IfI)=TYY<IfI)+BIG 
TYY<Ir J)«TYY<Ir J)-BIG 
TYY< Jr I )“TYY< JfD-BIG 
TYY< JfJ)=TYY<JfJ)+BIG 
SYY<IfI>=SYY(IfI)+BIG 
SYY<IfJ)“SYY<I*J)“BIG 
SYY<JrI)»SYY(JrI>~BIG 
JYY<Jr J)=SYY<Jr J)+BIG 
GO TO 500 
400 CONTINUE 
500 CONTINUE 
RETURN 
END 


SUBROUTINE STIFF for Example Problem 6 
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Example 7. Uppor Mantle Convection 

The final example illustrates convection in the upper mantle. 

Figure 12 shows the element layout and dimensions for the analysis. 

Because of the reduced depth of convection, the viscosity was reduced 
23 2 

to 0.3 x 10 N*S/m“ which gave a Rayleigh number equal to 
R = 3948. 

The analysis was not completed to convergence and after five iterations 
showed signs it might diverge. This could easily be corrected by the 
addition of more elements, particularly in the radial direction. 

Figure 12b, which illustrates the isotherms obtained on the fifth itera- 
tion, does indicate, however, what one would expect, that is, convection 
cells with aspect ratios approximately equal to one. It is interesting 
to note that the maximum surface velocity for this case was found to 
be 0,34 cm/year, which compares favorably with actual plate velocities. 
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CONCLUSION 

Tho provious examples domonstrato the use o£ MANTLE to solve many 
of tho complex questions related to mantle convection. In those 
OAamplos, tho coarsonoss of tho mesh adds a significant constraint on 
tho flow. Actual investigations into particular problems should bo 
made with a much finer mesh. Also, tho oxamplos did not demonstrate 
many of tho features of MANTLE, particularly transient problems 
including the use of quasi- Lagrangian analysis. The user is referred 
to the litoraturo cited in Appondix D for many more oxamplos covering 
a widor scope of application. 

As with all computer codes, improvements are continually being 
incorporated into MANTLE, At the timo of this roport, a now version 
of the code which includes Inertia terms in the equations c£ motion, 
the ability to model compressible elastic solids in the elasto- 
vi scop las tic analysis, an improved iterative technique for elasto- 
viscoplastic analysis, and improved solution algorithms for the matrix 
equations are being incorporated. These additions are not likely to 
significantly chango the input/output variables as described in this 
report. New versions of the code will be obtainable from the author 
when ready. 

gfiEGEDtWO F aGE BLANK N ° r PlL 
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APPENDIX A 

MESH2 GENERATOR AND IORDER ARRAY 
MGSH2 K«SW»» PAGE IW NOI HIM 

The basic building block of MESH2 is a quadrilateral region 
referred to as a loop. Any number of loops may bo joined together to 
form a completed mesh. Figure 1 illustrates how four such loops wore 
joined together to form the completed mesh shown. The four sides of a 
loop are parabolas specified by two corner points (dark circles in 
Figuro 2) and side point (open circles in Figuxe 2). 

The eight points specifying the shape of a given loop must begin 
with a comer point and be read in counterclockwise. The first three 
points then define side 1 of that loop with the other sados numborod 
consecutively counterclockwise. The user specifies how many divisions 
are desired on side 1 and side 2 of each loop (NDIV-array) . By use of 
the coordinates of the side points, the divisions (or elements) may 
then be "pushed" to one side or the other, e.g., compare the elements 
of loop 3 and the location of the side points for that loop. 

Loops are joined together through the use of array- JOIN. The data 
cards for JOIN are listed below for the four loops shown in Figure 1 
and Figure 2. Note, that each card only specifies how the present loop 
is joined to already existing loops. 
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Example of a Completed Mesh 



Exploded View Showing the 4 
Loops Used in MESH 


Figure 1 
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The two numbors listod for each side represent the loop and side 
respectively to which that side is joined. Hence, the above cards 
state for loop 1, thoro are no previous loops to which it is joined 
(this is always the case ] ; for loop 2, its sida 3 is joined to loop 1, 
side 1; for loop 3, its side 1 is joined to loop 2, side 2 and its 
side 4 is joined to loop 1, side 2; and for loop 4, its side 1 is 
joined to loop 1, side 3 and its side 2 is joined to loop 3, side 3. 

Figure 3 indicates the nodal point numbering and element numbering 
assigned by MESH for a given loop. For loops greater than the first, 
thoso numbers simply take up where the last previous loop left off. 

Where a side is joined to a previous loop, the previous loop numbers are 
used and the numbers of the current loop are adjusted so that they 
remain in consecutive order. 

The diagonals within each "square" are taken as the shortest 
distances. Wren those distances are equal, the direction of the 
diagonal is alternated as shown. 
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Figure 3 
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Nodal Point Numbering and tho IQKDER Array 

In order to keep tho bandwidth of tho stiffness matrix as small as 
possible it is possible to specif/ the order in which the elements will 
be selected. Some general guidelines are now given which will help the 
usor soloct an efficient ordering, 

At any stage of a frontal analysis it is always possible to draw 
one or more lines across the mesh which separates the elements not yet 
assembled from those already assembled. It is the number of nodal 
points along these lines which generally determine the bandwidth neces- 
sary for the analysis. It is therefore desirable to select an element 
order so that the maximum number of nodes at any stage in the frontal 
process is a minimum. Figure 4 shows a sequence of such lines, a-a 
through m-m, for the sample mesh. Although none of the lines shown 
happen to be the one where the maximum number of points will occur, 
they do give a good visual sensation of the "front 1 ' passing through 
the mesh. Lines d-d through lines j-j all have 19 nodal points along 
them. As tho front moves from the line g-g, so that it picks up the 
next element on the left between line g-g and h-h, four new nodal 
points will be added without a single loss. Hence, in this case the 
maximum number of points appears to be 23 rather than 19. The easiest 
way to select the path of the front is to locate what would appear to 
be the line passing through the mesh which would contain the largest 
number of points and yet be an obvious front for the analysis. In 
the example this was assumed to bo line g-g, From this line, it is a 
rather relatively easy matter to work outward from there so as to cover 
the entire mesh with a sequence of such lines. In the above example 
the I0RDER proceeded from the element at the beginning of row A, 
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n ntd t m 



Figure 4 
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through the four olomonts of that row and then on to the eight elements 
in row B beginning at tho loft and so forth through the mesh. It is 
important to emphasizo that this order does not influonco tho nodal 
point or element numbering gonoratod by MESH2. 
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APPENDIX B 

QUADRATURE POINTS, NODAL POINTS, AND STRESS POINTS 
IN THE 5-n PLANE 

The following figuro and table locato the quadrature points, 
nodal points, and stross points in tho £-n piano Chore, tho £ and 
n coordinates also correspond to two of the three area coordinates 
commonly referred to with Gaussian quadrature) . In order to locate a 
particular quadrature point or stress point on an element in the x-y 
plane, it is only necessary to orient tho nodal point numbers in the 
g-n plane with those in the x-y plane and assume the same relative 
orientation for the quadrature points and the stress points. 
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n, m 




n 



(b) Quadrature Points 


Quadrature 

Point 

5 

n 

Stress 

Point 

1 

0.33333 

0,33333 


2 

0.0S972 

0.47014 


3 

0.47014 

0.05972 


4 

0.47014 

0.47014 


5 

0.10129 

0.79743 

1 

6 

0.10129 

0,10129 

2 

7 

0.79743 

0.10129 

3 


Figure 1. Nodal Points, Quadrature Points, and Stress Points 
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APPENDIX C 
USEFUL REFERENCES 

Tho following references pertain directly to the material 
presented in this report. 

Thompson, E. G. and M. I. Hague, "A High Order Finite Element for 

Completely Incompressible Creeping Flow," Int. Journal Num. Moth . 
Engr ., Vol. 6, 315-321, 1973, 

Sato, A. and b, G. Thompson, "Finite Element Models for Creeping 
Convection," Journal Computational Physics , 

Thompson, E, G., "Average and Complete Incompressibility in the Finite 
Element Method," Int. Journal Nuro, Moth. Engr ., Vol, 9, 

925-932, 1975. 

Dawson, P. R., Finite Element Thormomechonicol Models for Metal Forming , 
Ph.D. Dissertation, Colorado State University, 1976. 

Crandall, S. H., Engineering Analysis , McGraw-Hill, 1956. 

Dawson, P. R, and E. G. Thompson, "Steady State Thermomechanical 
Finite Element Analysis of Elastoviscoplastic Metal Forming 
Processes," Numerical Modeling of Manufacturing Processes , ASME, 
PVP-PB-02S, 1977. 

Dawson, P. R, and E. G. Thompson, "Finite Element Analysis of Steady- 
State Elasto-Visco-Plastic Flow by tho Initial Stress-Rate Method," 
Int. Journal Num. Meth, Engr ., Vol. 12, 1978. 

Dawson, P. R. and E. G, Thompson, "Addendum to the Paper Finite Element 
Analysis of Steady-State Elasto-Visco-Plastic Flow by the Initial 
Stress-Rate Method," Int, Journal Num. Moth. Engr ., Vol, 12, 
pp. 382-383, 1978. 
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APPENDIX D 

INDEX OP FORTRAN VARIABLES 


Variable 

Pttf|0 

BNDRY 

69 

BNDRYC 

69 

BNDRYT 

69 

CII 

32,37,52 

CIJBC 

37 

GUI 

30 

CHO 

30 

COSLBC 


COSXXP 

32,37,5? 

CP 

68,76 

CP TIME 

57 

CTEMP 

46 

ex 

32,37,52 

CXI 

30 

CXLBC 

37 

CXO 

30 

CY 

32,37,52 

CYI 

30 

CYLBC 

37 

CYO 

3D 

DELP 

57,58 

DELQ 

58 

DELT 

57,58 

DELI! 

57,58 

DFCONV 

48 

DTCONV 

48 

DTIMB 

58 

DTMAX 

45 

DQCONV 

48 

DUCONV 

48 

DUMAX 

45 

EPSII 

67 

PTX 

59 

FTY 

59 

G 

67,76 

GAMX 

66,76 

GAMY 

66,76 

IB 

54 

IBMAX 

56 

IBMIN 

39 

ICOMP 

54 

IELE 

54 

IELEJ 

68,70 

IELEX 

34 


Variable 

PaRe 

IEMPT 

54 

I FLOW 

26 

INCPL 

42 

INCPR 

42 

INCPU 

42 

INCR 

50 

INTEMP 

44 

INTLCU 

47 

INTLTU 

47 

INTO 

54 

INTPL 

43 

INTPR 

43 

INTPU 

43 

^RDER 

40 

1‘4'NCII 

29,34 

pq 

60 

IREAD 

40 

/RZ 

44 

ITERC 

57,58 

ITERT 

57,58 

ITMAXC 

46 

ITMAXT 

46 

ITV 

44,70 

JBGN 

50 

JBND 

50 

JOIN 

35 

KMAX 

56 

KVOL 

54 

LAGEUL 

44 

LCU 

47 

LIST 

54 

LISTX 

54 

LPBC 

37 

LTD 

47 

MAT 

42 

MAXVOL 

39 

MJ 

66,67,68 

MOP 

44 

MOVE 

54 

MOVEX 

54 

MSHADJ 

69 

MSHCD 

39 

MNI 

45 

NDIV 

35 

NDIvR 

29 
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Variable 


V: "iablo 

PaRQ 

NDIVTil 

29 

TX 

32,37,52 

NCLMC 

51 

TXI 

30 

nelmt 

51 

TXO 

30 

NP 

S3 

TXLBC 

37 

NPBC 

32,37,52 

TY 

32,37,52 

NPBCI 

30 

1 i’lPC 

37 

NPBCO 

NPPG 

NPR 

NPTS 

NQMAX 

NSEC 

30 

27,33,67 

54 

49 

56 

49 

UT 

UX 

uxo 

UY 

UYO 

49,59 

49,59 

50 

49,59 

50 

NSDG 

54 

VECTL 

46 

NUMAT 

42 

vise 

67,76 

NUMBC 

31,36,37 

VS 

67,76 

NUMLPS 

33 

VT 

67,76 

NUMNP 

NUMPP 

NUMSEG 

NUMTP 

NUMVP 

51 

51 

56 

51 

51 

XBC 

XBC1 

XBCO 

XCOU 

XK 

32,37,52 
30 
30 
36 
66 , 67 

PENLTY 

67 

XLBC 

37 

PRS 

60 

XMAX 

34 

QRS 

60 

XMIN 

XORD 

34 

49,52,59 

RHO 

R1 

RKX 

RKY 

RM 

RPI 

RO 

68 

29 

68,76 

68,76 

29 

29 

20 

YBC 

YBCI 

YBCO 

YCOR 

YK 

YLBC 

YMAX 

32,37,52 

30 

30 

36 
67,76 

37 
34 

SIGH 

61 

YMIN 

34 

SIGXX 

61 

YORD 

49,52,59 

SIGXY 

61 



SIGYY 

61 



STIFF 

70,92 



TBC 

32,37,52 



TBCI 

30 



TBCD 

30 



TEMPK 

66,67 



TEMPO 

50 



THETA 

44 



TIMEM 

45 



TLBC 

37 



TRANS 

44 



TQ 

32,37,52 



TQI 

30 



TQO 

30 



TQLBC 

37 





